First-Principles Studies of the Metallization and the Equation of State of Solid Helium 
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The insulator-to-metal transition in solid helium at high pressure is studied with first-principles 
simulations. Diffusion quantum Monte Carlo (DMC) calculations predict that the band gap closes 
at a density of 21.3 g/cm 3 and a pressure of 25.7 terapascals, which is 20% higher in density and 
40% higher in pressure than predicted by density functional calculations based on the generalized 
gradient approximation (GGA). The metallization density derived from GW calculations is found 
to be in very close agreement with DMC predictions. The zero point motion of the nuclei had no 
effect on the metallization density within the accuracy of the calculation. Finally, fit functions for 
the equation of state are presented and the magnitude of the electronic correlation effects left out 
of the GGA approximation are discussed. 
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At low pressure, helium is an inert gas that exhibits 
a very large electronic excitation gap of 19.8 eV and has 
the highest ionization energy of all atoms, 24.6 eV. This 
is because helium has no core electrons, so its valence 
electrons are bound more strongly than in heavier atoms 
where screening effects play a role. Given such a strong 
binding, extreme pressures are needed to reach metalliza- 
tion. In fact, after neon pJQ], solid helium is expected 
to have the highest metallization pressure among all ele- 
mental solids. 

Metallic solid helium is expected to be present in the 
outer layers of white dwarfs (WD) 0]. After the initial 
star has exhausted all its nuclear fuel, it sheds its outer 
layer and leaves behind a dense carbon-oxygen core of 
the size of the earth that is surrounded by an envelope of 
pure helium, hydrogen, or a mixture. The fossil star then 
spends the remaining of its lifetime cooling until vanish- 
ing luminosity. Measuring the luminosity of the oldest 
WD would therefore constrain the age of the galaxy 
which qualifies WDs as stellar chronometers. 

Extracting the correct physics from WDs depends on 
how consistent the cooling models are with the observed 
luminosity [f| . Characterizing helium at high pressure is 
important because its properties regulate the heat trans- 
port across the outer layers. The metallization transition 
is important because it marks the point where the heat 
transfer switches from electronic conduction in interior 
WD layers to photon diffusion in the exterior. 

Most WD models rely on semi-analytical descriptions 
in the chemical picture Q where one treats helium as a 
collection of stable atoms, ions, and free electrons inter- 
acting via approximate pair potentials. While such ap- 
proaches work well at low density, they cannot describe 
the complex interactions in a very dense system, and a 
more fundamental description is required instead. First- 
principle methods, such as density functional theory and 
diffusion quantum Monte Carlo, are necessary as they 
provide a full accounting of quantum and statistical laws 
that govern the electrons and nuclei. 

Recent calculations by Kietzmann et al. [6j relied 
on DFT and the generalized gradient approximation 
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FIG. 1: Comparison of band gaps as function of density. DMC 
results for a small 2x1x1 rectangular supercell lie parallel and 
above the GGA gaps by 1.4 eV. DMC (4x2x2 rectangular cell) 
and GW metallization density of 21.3(1) g/cm 3 are in agree- 
ment. GGA underestimates the band gap by about 4 eV. The 
inset shows a part of the electronic band structure including 
the indirect gap. The filled and open circles indicate, respec- 
tively, the highest occupied valence state at the M point and 
the lowest unoccupied conduction state at F. Two densities 
are shown. The insulating state (black curve) lies below the 
metallic state (red dashed curve). 



(GGA) [7] to calculate the electrical conductivity and 
locate the insulator-to-metal transition in dense fluid he- 
lium. Kowalski et al. Q went beyond GGA to study the 
EOS, the electrical and optical properties of fluid helium 
up to 2 g/cm 3 . Stixrude and Jeanloz [§] studied the band 
gap closure in fluid helium over a wide range of densities 
including conditions of giant planet interiors. 

In this Letter, we study the metallization in solid 
helium at densities above 2 g/cm 3 using several first- 
principle simulation techniques. Our intention is to give 
an assessment of the accuracy of the widely used GGA 
for calculating total energies and band gaps at extreme 
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conditions. For this reason, we use the accurate, but ex- 
pensive diffusion quantum Monte Carlo (DMC) method 
and compare with the GW approximation to correct the 
GGA band gaps. We also derive the equation of state 
(EOS) and study the effect of the zero point motion of 
the nuclei on the band gap closure. We expect our first- 
principle study to serve as a guide for future laser experi- 
ments that are planned to extend the EOS measurements 
to very high pressures [IoT |. 

Our DFT calculations were performed with the 

ABINIT plane- wave basis code [llj. The electron- nuclei 
interactions were treated by a local Troullier-Martin 
norm-conserving pseudopotential 12] with a core radius 
of 0.4 a.u. We use Perdew-Burke-Ernzerhof GGA [7|] for 
the exchange-correlation functional. We worked with an 
8x8x8 Monkhorst-Pack k-point grid and a plane wave en- 
ergy cutoff of 230 Ha. 

Mao et al. have demonstrated experimentally that 
solid helium still remains in the hep structure up to high 
pressures (57 GPa), apart from a limited fee loop along 
the melting line at low temperatures. So we adopt the 
hep solid phase and optimize the cell geometry iteratively 
until the pressure has converged to within 1%. 

The band structure in the inset in Fig. Q] shows that 
solid helium in the hep structure exhibits an indirect 
band gap. The lowest unoccupied state occurs at the 
r wave vector. The highest occupied state is at k = 
0.95kM, which is very close to the M wave vector. Sub- 
sequently, we approximate the excitation gap by the dif- 
ference in energy between the valence M point and the 
conduction T point. As Fig. [T] shows, the band gap de- 
creases almost linearly with density. GGA predicts the 
band gap closure at a density of 17.4 g/cm 3 , which cor- 
responds to a pressure of 17.0 TPa. 

Kohn-Sham GGA is known to systematically underes- 
timate the band gap. DMC is expected to predict the 
band gap width more accurately, because it explicitly in- 
cludes electronic correlation effects. In fact, DMC has 
been used successfully to describe the electronic ground 
state in weakly as well as in strongly correlated sys- 
tems Excited states were also calculated reliably 
with DMC @, One needs a large supercell to describe 
all correlation effects, which comes at a high computa- 
tional cost since DMC scales as high as 0(aN 3 + bN 2 ), 
where N is the number of electrons. 

Our DMC simulations were performed with the 
CASINO code [17]. In DMC, the Schrodinger equation 
is solved stochastically by simulating branching and dif- 
fusion in imaginary time. A trial wave function enters in 
the propagation of electronic configurations and we use 
the fixed-node approximation to avoid the fermion sign 
problem. Our trial wave function is of the Slater-Jastrow 
form. The parameters in the Jastrow factor were opti- 
mized by variance minimization. They comprise electron- 
electron, electron-nucleus and electron-electron-nucleus 
terms, that the We use GGA orbitals for the Slater part. 



We also keep the same pseudopotential as in DFT calcu- 
lations. We picked a conservative high energy cutoff of 
800 Ha. For efficiency, the orbitals are represented nu- 
merically using blip functions [iH ]. Our results are well 
converged with a time step of 0.002 a.u. To calculate 
the band gap, we promote one electron, either spin-up 
or -down, from the valence band to the conduction band 
(M — ► T). The calculations for the excited state used 
the same Jastrow parameters as for the ground state cal- 
culation because the DMC energy depends only on the 
nodes of the trial wavefunction that are determined by 
the Slater determinant. To include the Y and M wave vec- 
tors into one DMC band gap calculation, we have chosen 
rectangular 2x1x1 or 4x2x2 supercells with, respectively 
16 or 128 electrons. We also performed simulations with 
up to 3x3x3 triangular supercells with 108 electrons to 
determine the ground state energies and to study their 
finite-size dependence. 

Figure [T] shows the DMC band gaps to be larger than 
the GGA results, as expected. The DMC curves for 
the 2x1x1 and 4x2x2 rectangular supercells lie parallel 
to the GGA curve, hence showing a gap correction that 
is independent of density. DMC simulations in a 4x2x2 
rectangular cell predict a metallization density of 21.3(1) 
g/cm 3 . We consider the 4x2x2 gap results to be con- 
verged with respect to system size because the ground 
state energy agrees well with triangular 3x3x3 supercells 
and the remaining finite size corrections are small. 

The 4x2x2 DMC band gaps agree very well with the 
GW results, which also show a linear behavior with den- 
sity similar to all previous calculations. The GW approx- 
imation has proven to be a reliable method for correcting 
the GGA band gaps in a variety of materials [l9j . Our 
GW band gap corrections are calculated within an accu- 
racy of 0.1 eV after converging the number of bands (50) 
and plane waves (27 Ha). In comparison, our metalliza- 
tion density is significantly higher than the linear-muffin- 
tin-orbitals prediction of 13.5 g/cm 3 for helium in the fee 
phase (20j . 

We tried improving the nodes of the trial wave function 
in DMC by adding a backflow correction to the Slater de- 
terminant. This method introduces further correlation to 
the trial wave function by replacing the electron coordi- 
nates in the determinant by a set of collective coordi- 
nates. The DMC total energy decreased slightly but the 
band gap did not change within error bars (0.02 eV). 

We also studied whether the zero point motion of the 
nuclei has any effect on the metallization density because 
one might expect the disorder introduced by the zero 
point motion to reduce the metallization density, as al- 
ready noted in the fluid 8, 9] . We generated series of con- 
figurations with path integral Monte Carlo (PIMC) [2l[ 
simulations for different temperatures between 500 and 
5000 K. The helium atoms interact with an effective pair 
potential that we constructed by matching the forces [28| 
of a density functional molecular dynamics (DFT-MD) 
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FIG. 2: Panel (a): Difference between DMC and GGA ground 
state energies. We calculate the DMC energy in two ways, 
by extrapolations to infinite cell size (circles) and by di- 
rectly using results from our largest 3x3x3 triangular su- 
percell (diamonds). In panel (b), we relate the DFT-DMC 
energy difference to the amount of mechanical work needed 
to reach a certain density. We plot the correlation fraction, 
(£ D FT-BDMc)/(BDFT-B a tom), where E atom = -79.0048 eV 
is the exact energy of the isolate helium atom. 



simulation at 5000 K. We verified the accuracy of the 
potential by comparing the original DFT-MD pair corre- 
lation function with that of classical Monte Carlo simu- 
lations. We then computed the GGA band gaps for the 
PIMC configurations and compare with perfect hep lat- 
tice results. Within error bars, the zero point motion had 
no effect on the band gap. 

DMC is computationally intensive which limits the size 
of the supercell used to simulate the infinite solid. Our 
biggest supercells consist of 3x3x3 triangular (108 elec- 
trons) and 4x2x2 rectangular (128 electrons) primitive 
unit cells. We correct the energies by considering the fol- 
lowing finite size effects. To first order, the independent- 
particle finite size effects (IPFSE) dominate [22|]. The 
error arises from an incomplete k-point sampling of the 
Brillouin zone in the DMC supercell. In DFT, the com- 
putational cost is directly proportional to the number of 
k-points because the electrons are treated as indepen- 
dent particles. In DMC, however, bigger supercells are 
needed to include more k-points. So the limitation is 
more severe since the computational cost in DMC scales 
with the number of particles in the supercell as 0(aN 3 + 
bN 2 ). Our IPFSEs range from 0.07 eV/el up to 2.0 eV/el 
with increasing density in our 3x3x3 triangular cell. We 
also include kinetic energy corrections of the order of 
0.1 eV/el, which are due to long range correlations [23| . 

We correct for Coulomb finite size effects (CFSE) that 
are due to the long range interaction between charged 
particles and their periodic images. These effects decay 
with system size as 1/N and tend to lower the energy 
slightly. We reduce these effects by using the model in- 
teraction potential (MPC) instead of the Ewald interac- 
tion The difference between the total ground state 



energies computed with these two interactions is 0.04 
eV/el for the largest volume. The error increases with 
density because the periodic image charges are closer in 
smaller supercells. We obtain very similar band gaps 
with the Ewald and the MPC interactions. This is due 
to the cancellation of CFSE errors when taking energy 
differences to calculate the band gap. 

After correcting the DMC ground state energies for 
finite size effects, we compare with DFT calculations in 
Fig. [2] We use two ways to estimate the correction energy 
that is missing in GGA. First we report the DFT-DMC 
difference directly for our larger 3x3x3 supercell. Sec- 
ondly we extrapolate the DMC energies to infinite size 
as function of 1/N. We derive an uncertainty of the re- 
sulting DMC energies by comparing linear and quadratic 
extrapolation. The resulting error bars in Fig.[5]are small 
except for the highest density. In the density range of 
consideration, the correlation energy error in GGA is ap- 
proximately constant, 0.36 eV/el. 

The correlation error in GGA becomes less important 
with increasing density because it stays constant while 
the energy increases with compression. This trend is il- 
lustrated in Fig. [5] where we relate the DFT-DMC energy 
difference to the energy needed to compress helium to 
density, p. The plotted ratio, (-Bdft — -Edmc)/(-£'DFT — 
-E"atom)) approaches unity in the low density limit. In 
the opposite high density limit, the graph shows how 
it tends to zero as helium approaches the state of a 
one-component plasma with a neutralizing background. 
The kinetic energy of homogeneous electron gas domi- 
nates the correlation and eventually the Coulombic en- 
ergy terms. GGA is expected to describe this limit 
well since the exchange-correlation functional was de- 
rived from DMC simulations of the homogeneous electron 
gas 0. 

Since the corrections to the GGA energy appear to 
be independent of density, there are no corrections to 
pressures derived from GGA. We were able to represent 
our zero-temperature static lattice GGA pressure data 
in the insulating regime by a Vinet EOS curve 25[ with 
the parameters Vq = 20.6397 cm 3 /mol as zero-pressure 
volume, B = 0.01928 GPa as bulk modulus, and B' Q = 
9.2153 as its derivative. The fit reproduced our GGA 
data from 3 to 1200 GPa with an accuracy of 3%. The 
comparison with experiments [2(| was studied in great 
detail earlier [27 1. 

It was not possible to extend the Vinet fit into the 
metallic regime. Instead, we adopt a fit based on the 
parametrization of the homogeneous electron gas energy. 
This fit includes the kinetic, Coulombic, exchange as well 
as correlation terms and, in this case also ionic contribu- 
tions, P(V) = ^3+^3+^ + ^3. In units of 
GPa and cm 3 /mol, the coefficients are a\ = 3186.21, 
a 2 = -2761.74, a 3 = -565.78, and a 4 = 854.71 where 
the leading coefficient is taken from the free Fermi gas. 
The fit reproduces our DFT data points from 1.2 to 1600 
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FIG. 3: Pressure-volume relationship derived from GGA cal- 
culations for solid helium. The experimental Vinet fit [2rj| . 
valid up to 57 GPa, agrees with the low pressure Vinet fit 
of GGA data. The high pressure (HP) fit to the GGA data 
approaches the free homogeneous electron gas (HEG) . To em- 
phasize differences, PV 3 was plotted in the main graph. 



TPa within 0.5%. In Fig. [3l we show the pressure over a 
large density range. In the high density limit, the corre- 
lation effects decrease and the DFT pressure approaches 
the free homogeneous electron gas behavior. 

In conclusion, we have demonstrated that solid helium 
reaches a metallic state at an extreme pressure of 25.7 
TPa, which is significantly larger than predicted by stan- 
dard GGA method. For WD interiors, this implies that 
the inner layer of metallic helium is thinner and the outer 
region where photon diffusion dominates the heat trans- 
port is larger than previously predicted. 

With quantum Monte Carlo we have shown that stan- 
dard GGA methods underestimate the band gap in solid 
helium by 4 eV, which translates into an underestimation 
of the metallization pressure by 40%. The GW band gap 
corrections are in good agreement with DMC calcula- 
tions, which offers the possibility of using GW for cor- 
recting the band gaps derived from GGA simulation of 
fluid helium at high pressure and to make more realistic 
comparisons with shock wave measurements of conduc- 
tivity and reflectivity. 

Finally, we determined the equation of state and pre- 
sented a fit. We analyzed the correlation effects that are 
missing in GGA and demonstrated that their importance 
decreases relative to the total energy with increasing den- 
sity as helium approaches the state of a one-component 
plasma with a rigid neutralizing background. 

This work was supported by NSF and NASA. We 
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